%% Description
% In this introductory example, a block of "standard tissue" (mu_a = 1,
% mu_s = 100, g = 0.9) is illuminated by a pencil beam (infinitely thin
% beam). A small slice of air is present in the top of the simulation
% volume. nx and ny are set to odd values so that when the pencil beam is
% launched at x = y = 0 and travels straight down, it travels along a
% well-defined center column of voxels (the middle of the 51st column). Use
% the log10 plot checkbox in the visualizations to better see the fluence
% rate and absorption distribution in the MC result.

%% MCmatlab abbreviations
% G: Geometry, MC: Monte Carlo, FMC: Fluorescence Monte Carlo, HS: Heat
% simulation, M: Media array, FR: Fluence rate, FD: Fractional damage.
%
% There are also some optional abbreviations you can use when referencing
% object/variable names: LS = lightSource, LC = lightCollector, FPID =
% focalPlaneIntensityDistribution, AID = angularIntensityDistribution, NI =
% normalizedIrradiance, NFR = normalizedFluenceRate.
%
% For example, "model.MC.LS.FPID.radialDistr" is the same as
% "model.MC.lightSource.focalPlaneIntensityDistribution.radialDistr"

%% Geometry definition
MCmatlab.closeMCmatlabFigures();
model = MCmatlab.model;

model.G.nx                = 201; % Number of bins in the x direction
model.G.ny                = 201; % Number of bins in the y direction
model.G.nz                = 40; % Number of bins in the z direction
model.G.Lx                = 2; % [cm] x size of simulation cuboid
model.G.Ly                = 2; % [cm] y size of simulation cuboid
model.G.Lz                = .4; % [cm] z size of simulation cuboid

model.G.mediaPropertiesFunc = @mediaPropertiesFunc; % Media properties defined as a function at the end of this file
model.G.geomFunc          = @geometryDefinition; % Function to use for defining the distribution of media in the cuboid. Defined at the end of this m file.

model = plot(model,'G');

%% Monte Carlo simulation
model.MC.simulationTimeRequested  = .1; % [min] Time duration of the simulation
model.MC.matchedInterfaces        = false; % If false, uses the refractive indices as defined in mediaPropertiesFunc at the end of this file
model.MC.boundaryType             = 1; % 0: No escaping boundaries, 1: All cuboid boundaries are escaping, 2: Top cuboid boundary only is escaping, 3: Top and bottom boundaries are escaping, while the side boundaries are cyclic
model.MC.wavelength               = 808; % [nm] Excitation wavelength, used for determination of optical properties for excitation light

model.MC.lightSource.sourceType   = 0; % 0: Pencil beam, 1: Isotropically emitting line or point source, 2: Infinite plane wave, 3: Laguerre-Gaussian LG01 beam, 4: Radial-factorizable beam (e.g., a Gaussian beam), 5: X/Y factorizable beam (e.g., a rectangular LED emitter)

% For a pencil beam, the "focus" is just a point that the beam goes
% through, here set to be the center of the cuboid:
model.MC.lightSource.xFocus       = 0; % [cm] x position of focus
model.MC.lightSource.yFocus       = 0; % [cm] y position of focus
model.MC.lightSource.zFocus       = model.G.Lz/2; % [cm] z position of focus

model.MC.lightSource.theta        = 0; % [rad] Polar angle of beam center axis
model.MC.lightSource.phi          = 0; % [rad] Azimuthal angle of beam center axis


% These lines will run the Monte Carlo simulation with the provided
% parameters and subsequently plot the results:
model = runMonteCarlo(model);
model = plot(model,'MC');
% 提取结果数据
layer = squeeze(model.MC.normalizedFluenceRate(:,:,35));
%layer_log = log(layer);
figure;
imagesc(layer);
colorbar; % 添加颜色条
axis equal; % 设置坐标轴比例一致


% 傅里叶变换
figure;
layer_fft = fftshift(fft2(layer));
layer_fft_abs = abs(layer_fft);
imagesc(layer_fft_abs);
colorbar;
axis equal;
title('频谱图');
layer_fft_abs_avg = layer_fft_abs;
% 求傅里叶变换的圆心
[max_value, linear_index] = max(layer_fft(:));
[x, y] = ind2sub(size(layer), linear_index);
[row, col] = size(layer_fft);
distances_squared  = sqrt(((1:row).' - x).^2 + ((1:col) - y).^2);
distances = round(sqrt(distances_squared), 2);
%disp(distances);
unique_distances = unique(distances(:));
% 对每个距离值进行循环处理
for i = 1:numel(unique_distances)
    % 获取当前距离值
    current_distance = unique_distances(i);
    % 在 distances 中找到与当前距离值相等的位置
    positions = find(distances == current_distance);
    % 从 layer_fft 中获取对应位置的值，并计算均值
    mean_value = mean(layer_fft_abs(positions));
    % 将均值赋值给对应位置的 layer_fft
    layer_fft_abs_avg(positions) = mean_value;
end
figure;
imagesc(layer_fft_abs_avg);
colorbar;
axis equal;
title('频谱图平滑');
row_data_layer_fft_abs = layer_fft_abs(x, :);
row_data_layer_fft_abs_avg = layer_fft_abs_avg(x, :);
figure;
plot(row_data_layer_fft_abs);
hold on; 
plot(row_data_layer_fft_abs_avg);
hold off; 
% 设置图例
legend('layer\_fft\_abs', 'layer\_fft\_abs\_avg');
% 设置标签和标题
xlabel('position');
ylabel('amplitude');
title('FFT curve');
filename = 'mua0@1_mus20@8_g0@59.mat';
save(filename, 'layer_fft_abs', 'layer_fft_abs_avg', 'row_data_layer_fft_abs', 'row_data_layer_fft_abs_avg');
%% Geometry function(s) (see readme for details)
% A geometry function takes as input X,Y,Z matrices as returned by the
% "ndgrid" MATLAB function as well as any parameters the user may have
% provided in the definition of Ginput. It returns the media matrix M,
% containing numerical values indicating the media type (as defined in
% mediaPropertiesFunc) at each voxel location.
function M = geometryDefinition(X,Y,Z,parameters)
  zSurface1 = 0.05;
  zSurface2 = zSurface1 + 0.3; % 0.3 cm in thickness
  M = ones(size(X)); % Air
  M(Z > zSurface1) = 2; % "Standard" tissue
  M(Z > zSurface2) = 1; % Air
end

%% Media Properties function (see readme for details)
% The media properties function defines all the optical and thermal
% properties of the media involved by filling out and returning a
% "mediaProperties" array of "mediumProperties" objects with various
% properties. The j indices are the numbers that are referred to in the
% geometry function (in this case, 1 for "air" and 2 for "standard tissue")
% See the readme file or the examples for a list of properties you may
% specify. Most properties may be specified as a numeric constant or as
% function handles.
% 
% The function must take one input; the cell array containing any special
% parameters you might specify above in the model file, for example
% parameters that you might loop over in a for loop. In most simulations
% this "parameters" cell array is empty. Dependence on wavelength is shown
% in examples 4 and 23. Dependence on excitation fluence rate FR,
% temperature T or fractional heat damage FD can be specified as in
% examples 12-15.
function mediaProperties = mediaPropertiesFunc(parameters)
  % Always leave the following line in place to initialize the
  % mediaProperties array:
  mediaProperties = MCmatlab.mediumProperties;

  % Put in your own media property definitions below:
  j=1;
  mediaProperties(j).name  = 'air';
  mediaProperties(j).mua   = 1e-8; % Absorption coefficient [cm^-1]
  mediaProperties(j).mus   = 1e-8; % Scattering coefficient [cm^-1]
  mediaProperties(j).g     = 1; % Henyey-Greenstein scattering anisotropy
  mediaProperties(j).n     = 1;
  
  j=2;
  mediaProperties(j).name  = 'standard tissue';
  mediaProperties(j).mua   = 0.1; % Absorption coefficient [cm^-1]
  mediaProperties(j).mus   = 20.8; % Scattering coefficient [cm^-1]
  mediaProperties(j).g     = 0.59; % Henyey-Greenstein scattering anisotropy
  mediaProperties(j).n     = 1.49; %PMMA at 808 nm
  
  j=3;
  mediaProperties(j).name  = 'air';
  mediaProperties(j).mua   = 1e-8; % Absorption coefficient [cm^-1]
  mediaProperties(j).mus   = 1e-8; % Scattering coefficient [cm^-1]
  mediaProperties(j).g     = 1; % Henyey-Greenstein scattering anisotropy
  mediaProperties(j).n     = 1;
  
end